library(stringr)
library(car)
library(plyr)
library(survey)
library(ltm)
library(interplot)
setwd("~/Application3/")

## Configuring the NFC adaptive battery ----------------------------------------
cognition_total <- read.csv("data/NFCTrainingData.csv", row.names = 1)
nfc_model <- grm(cognition_total[ ,1:18], Hessian = TRUE, control = list(GHk = 21))

## Results in manuscript summarizing AI -------------------------------------------------------------
predictions <- factor.scores(nfc_model, method = "EAP",
                             resp.patterns = cognition_total[ ,1:18])

## 'sample' is turk1, turk2, or taps; 'score' is prediction
sample_index <- c(rep(0, 1204), rep(1, 1333), rep(2, 1506))

## reverse coding theta estimates so higher values mean "more" NFC
## (first discrimination param in 'nfc_model' is constrained to be positive
## and in this context that would mean that lower theta values indicate
## "more" NFC) 
pred_mat <- cbind(sample_index, -1*predictions$score.dat$z1)
colnames(pred_mat) <- c("sample", "score")
pred_mat <- data.frame(pred_mat)


## reading in tree created for ANES
nfc_tree <- read.csv("data/nfcTable.csv", stringsAsFactors = FALSE, row.names = 1)

## reading in pilot results
all_pilot_data <- read.csv("data/anes_pilot_2016.csv")

## Limit down to NFC items that are used in tree
## and recode missing to NA
nfc_items <- c("NFC1", "NFC10", "NFC15", "NFC21", "NFC23",
              "NFC24", "NFC29", "NFC32", "NFC40")
nfc_pilot_data <- all_pilot_data[ , nfc_items]
nfc_pilot_data[nfc_pilot_data == 9] <- NA

## Response rates for this sample... everyone answered 4
table(apply(nfc_pilot_data, 1, function(row) sum(!is.na(row))))


## number of respondents answering each item
num_answered <- apply(nfc_pilot_data, 2, function(c) sum(!is.na(c)))
sort(num_answered)

## number of complete branches including each item
complete_branches <- nfc_tree[apply(nfc_tree[,1:18], 1, function(r) sum(!is.na(r)) == 3), ]
num_branches <- (apply(complete_branches[,nfc_items[-c(5)]], 2, function(c) sum(!is.na(c))) +
                table(complete_branches[,"NextItem"]))
sort(num_branches)

## function to check branching schemes used
check_use <- function(branch, respondent){
  test1 <- sum(nfc_tree[branch, nfc_items] ==
                 nfc_pilot_data[respondent, ], na.rm = TRUE)==3
  test2 <- !is.na(nfc_pilot_data[respondent, nfc_tree[,"NextItem"][branch]])
  output <- FALSE
  if(all(c(test1, test2))){
    output <- TRUE
  }
  return(output)
}

branch_used <- rep(NA, 1200)
for(i in 1:1200){
  results <- aaply(1:259, 1, check_use, respondent = i)
  if(any(results)){
    branch_used[i] <- which(results)
  }else{
    branch_used[i] <- NA
  }
}

## number of branching schemes used
length(table(branch_used))

## a check that each respondent used a possible scheme
sum(table(branch_used))


## NFC scores for ANES Sample --------------------------------------------------
## Cat object
library(catSurv)
nfc_cat <- grmCat(nfc_model)
nfc_cat@priorParams <- c(mean(pred_mat$score[pred_mat$sample == 2]), 1.2)
nfc_cat@estimation <- "EAP"
nfc_cat@selection <- "EPV"

theta_ests <- function(x){
  empty_answers <- nfc_cat@answers
  names(empty_answers) <- names(nfc_cat@difficulty)
  answers <- as.numeric(nfc_pilot_data[x, nfc_items])
  ## index of 9 of 18 items actually used
  index_vec=c(1, 3, 5, 8, 10, 11, 12, 14, 17)
  if(!all(nfc_items == names(nfc_cat@difficulty)[index_vec])){
    print("not set up right")
    stop
  }
  new_answers <- empty_answers
  new_answers[index_vec] <- answers
  nfc_cat@answers <- new_answers
  return(estimateTheta(nfc_cat))
}


## reverse coding theta estimates so higher values mean "more" NFC
all_pilot_data$nfc_estimates <- laply(1:1200, theta_ests) * -1

## descriptive statistics: survey weighted mean and standard deviation
survey_stats <- svydesign(id = ~caseid, data = all_pilot_data, weight = ~weight)
svymean(~nfc_estimates, design = survey_stats)
sqrt(svyvar(~nfc_estimates, design = survey_stats))

## TAPS sample
mean(pred_mat$score[pred_mat$sample == 2])
sd(pred_mat$score[pred_mat$sample == 2])


## Syrian refugee framing experiment -------------------------------------------
all_pilot_data$syrians_a[all_pilot_data$syrians_a == 9] <- NA
all_pilot_data$syrians_b[all_pilot_data$syrians_b == 9] <- NA
all_pilot_data$syrians_b[all_pilot_data$syrians_b == 8] <- NA
all_pilot_data$refugee <- rowMeans(cbind(all_pilot_data$syrians_a,
                                         all_pilot_data$syrians_b),
                                   na.rm = TRUE)

## Recoding some variables
all_pilot_data$pid7 <- recode(all_pilot_data$pid7, "8=NA")
all_pilot_data$ftmuslim <- recode(all_pilot_data$ftmuslim, "998=NA")
all_pilot_data$ideo5 <- recode(all_pilot_data$ideo5, "6=NA")
all_pilot_data$ideo5 <- recode(all_pilot_data$ideo5, "6=NA")
all_pilot_data$isis_troops <- recode(all_pilot_data$isis_troops, "8=NA")
all_pilot_data$rr1 <- recode(all_pilot_data$rr1, "8=NA")
all_pilot_data$rr2 <- recode(all_pilot_data$rr2, "8=NA")
all_pilot_data$rr3 <- recode(all_pilot_data$rr3, "8=NA")
all_pilot_data$rr4 <- recode(all_pilot_data$rr4, "8=NA")
all_pilot_data$rr <- with(all_pilot_data, rr1 - rr2 -rr3 + rr4)
all_pilot_data$gender <- recode(all_pilot_data$gender, "2=1; 1=0")
all_pilot_data$RAND_ISIS <- recode(all_pilot_data$RAND_ISIS, "2=1; 1=0")

## Models
cw_mod <- lm(refugee ~ RAND_ISIS + nfc_estimates + educ + pid7 + 
                  ftmuslim + ideo5 + gender + isis_troops + 
                  I(race==1) + I(race==2) + I(race==3) + rr,
                data = all_pilot_data, weights = weight)

cw_nfc_mod <- lm(refugee ~ RAND_ISIS*nfc_estimates + educ + pid7 + 
                ftmuslim + ideo5 + gender + isis_troops + 
                I(race==1) + I(race==2) + I(race==3) + rr,
                data = all_pilot_data, weights = weight)

## Table 5 results
summary(cw_mod)
summary(cw_nfc_mod)


## Figure 5: Interaction plot
library(interplot)
pdf(width = 5, height = 4, file = "results/Figure5.pdf")
interplot(m = cw_nfc_mod, var1 = "RAND_ISIS", var2 = "nfc_estimates") +
  # Add labels for X and Y axes
  xlab("Need for cognition (AI version)") +
  ylab("Effect of civil war frame") +
  # Change the background
  theme_bw() +
  # Add the title
  # ggtitle("Estimated effect of civil war frame \non opposition to syrian refugees") +
  theme(plot.title = element_text(face="bold")) +
  # Add a horizontal line at y = 0
  geom_hline(yintercept = 0, linetype = "dashed")
dev.off()
